Localised plumes in three-dimensional compressible 

magnetoconvection 
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Abstract 

Within the umbrae of sunspots, convection is generally inhibited by the presence of strong vertical 
magnetic fields. However, convection is not completely suppressed in these regions: bright features, known 
as umbral dots, are probably associated with weak, isolated convective plumes. Motivated by observations 
of umbral dots, we carry out numerical simulations of three-dimensional, compressible magnetoconvection. 
By following solution branches into the subcritical parameter regime (a region of parameter space in which 
the static solution is linearly stable to convective perturbations), we find that it is possible to generate a 
solution which is characterised by a single, isolated convective plume. This solution is analogous to the 
steady magnetohydrodynamic convectons that have previously been found in two-dimensional calculations. 
These results can be related, in a qualitative sense, to observations of umbral dots. 

1 Introduction 

Using modern instruments, such as the Solar Optical Telescope on board Hinode and the 1-metre Swedish 
Solar Telescope on La Palma, it is possible to make detailed observations of magnetic fields and convection at 
the surface of the Sun. Sunspots are the most prominent magnetic features on the solar surface. A typical 
sunspot consists of a central umbral region, surrounded by a complex filamentary penumbra. Umbral regions 
appear dark because their surface temperatures are ( typically) only 70 — 85% of the mean surface temperature 
of the non-magnetic photosphere (see, for example. iThomas fc Weiss . 2008). This reduction in temperature 
is due to the fact that the convective transport of heat is impeded within sunspot umbrae by the presence of 
strong, near vertical magnetic fields (which can often exceed 3000G) . 

Detailed observations of sunspot umbrae have shown that they are not uniformly da rk. In almost all 
sunspots, bright point-like structures can be observed 



these are known as umbral dots (jDanielsonl . I1964D . 
These bright features are warmer than thei r immediate surroundings, but are (gene r ally) cooler than the 
surrounding photosphere (see, for example, Sobotka fc Hanslmeier . 20051 : Kitai et all . 12007 ). It is difficult 
to determine the characteristic size of an umbral dot, althou gh these features are always small compared to 
the umbral diameter. In a recent study, iKitai et al. ( 2007 ) found that the umbral dots in one particular 
sunspot had typical diameters of approximately 220 — 350km, although a significant number appeared to 
be much smaller than this (possibly below the resoluti o n lim it for the Solar Optical Telescope on Hinode) . 
Umbral dots are also short-lived features. IKitai et al.l ( 2007 ) found that most of the umbral dots in their 



survey had lifetimes of between 5 and 20 minutes. In an earlier study, Sobotka et al. ( 1997t ) found a much 
bro ader range of lifet imes for umbral dots (with a small percentage lasting longer than 2 hours), although, 
like Kitai et al. | (|2007l) . they found a mean lifetime of approximately 15 minutes. Most umbral dots exhibit no 
systematic proper motions. However, those that appear to form at the umbral/penumbral boundary (which 
are often associated with penumbral gr ains) tend to migrate radially inwards towards the centre of the umbra 
( Sobotka et"ailll995HKitai et all 120071). There is som e observational evidence for weak upflows within umbral 



dots 



2007 



Socas-Navarro et all 2004 : iBharti et al. . 2007 ) as well as downflows around their edges ( Bharti et al 



Ortiz et alll201fj) . Clearly, the observations indicate that umbral dots correspond to convective plumes 
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"within sunspot umbrae. Further theoretical support for this conclusion comes from the work of lDeinzerl (|1965I) . 
who determined that convective motions must be present within the umbra, as radiative processes alone could 
not transport sufficient energy to the surface. 

Theoretical studies of umbral convection tend to be based upon local models of magnetoconvection in a 
Cartesian domain. It is well-kn own that a strong vert ical magnetic field tends to inhibit convective motions in 
an electrically-conducting fluid ([Chandrasekhaii Il96ll ). When the dynamics are dominated by ma gnetic fields 



conve ction takes the form of weak, narrow plumes. In an idealised model of magnetoconvection, Weiss et al 



(2002) bund a steady, almo s t hexa gonal pat tern of convection in the magnetically-dominated regime. More 
recently, ISchtissle~ fc Voglerl (l2006h (see also iBharti et all l201dh have carried out a more realistic set of cal- 
culations, including the effects of partial ionisation and radiative transfer. These simulations produced a 
time-dependent pattern of individual convective p lumes, the p r opert ies of which compare very favo urably to 
observations of umbral dots. In the calculations of Weiss et al. ( 2002 ) and Schussler fc Vogierl ( 20061 ). convec- 
tive features tend to be distributed across the whole computational domain. It is worth noting that observation s 
indicate that the distribution of umbral dots tends to be rather non-uniform ( Sobotka fc HanslmeierL 20051) . 
This could simply be explained by variations in intensity (as seen in the calculations of ISchiissler fc Vogler , 
2006), but some degree of localisation in the distribution of convective plumes could also help to explain the 



observed distribution of umbral dots. 

In a two-dimensional model of incom pressibl e mag netoconvection, iBlanchflower (1999) found strongly lo- 
calised, steady convective states (see also Dawesl . 120071) . These were named convectons. Restricting attention 
to a simplified model, in which the governing equat ions were projected onto a minimal set of Fourier modes in 
the vertical direction, IBlanchflower fc Weis s (2002) were also able to find oscillatory localised states in three 
spatial dimensions. In both these cases, these localised states were found in the subcritical parameter regime 
(in which the static, purely conducting state is linearly stable to convective pe rturbation s ) . Th ese localised 



states are an extreme example of a p henomenon th at is known as flux separation ( Tao et al. . 1998). Convective 



plumes tend to expel magnetic flux ( Weisd . 19661 ). causing it to accumulate in the surrounding fluid. In this 
subcritical parameter regime, the magnetic field that surrounds the plume becomes sufficiently strong that 
convection in the surrounding fluid is completely inhibited, giving rise to a truly localised convective state. If 
conditions are appropriate for subcritical convection within sunspots (something that is certainly plausible), 
these results suggest that it may be possible for truly localised umbral dots to form within sunspot umbrae. 

In this paper, we demonstrate the e xistence of steady loca l ised c onvective plumes in three-dimensional 
compressible magnetoconvection. Unlike IBlanchflower fc Weissl (J2002J), we make no simplifying assumptions 
regarding the vertical structure of the convective flows. The setup of the model is described in detail in the 
next section of the paper. Numerical results from this model are presented in Section 3. In the final section, 
we relate our findings (in qualitative terms) to observations of umbral dots. 



2 Problem description and setup 

We consider the evolution of a layer of compressible, electrically-conducting fluid, heated from below, in the 
presence of an imposed magnetic field. Various properties of the fluid, including the thermal conductivity, K, 
the shear viscosity, /x, the magnetic diffusivity, 77, the magnetic permeability, hq, and the specific heat capacities 
at constant pressure and density (cp and cy respectively) are assumed to be constant. At a position x and time 
t, we define p(x, t), T(x, i) and u(x, t) to be the fluid density, temperature and velocity field (respectively), 
whilst B(x, t) represents the magnetic field. 

This fluid occupies a three-dimensional Cartesian domain with < z < d and < x, y < 8d. The axes of 
this coordinate system are orientated so that the z-axis points vertically downwards, parallel to the constant 
gravitational acceleration, g = gz. For this model problem, periodic boundary conditions are imposed in the x 
and y directions, whilst the upper and lower boundaries (at z = and z = d) are assumed to be impermeable 
and stress free. Furthermore, fixed temperature boundary conditions are applied at the upper and lower 
boundaries with T = Tq at z = and T = Tq + AT at z = d (AT > 0). It is also assumed that the horizontal 
components of any magnetic fields that are present vanish at z — and z — d. When the layer is static, the 
imposed magnetic field is uniform and vertical, i.e. B = BqZ. 

Before writing down the governing equations fo r this system, w e can expr ess these in non-dimensiona l 
form. More details of this procedure can be found in lMatthews et aL ( 1995 ) and Bushbv fc Houghton! ( 20051 ). 



1/2 

Very briefly, all lengths are scaled by the layer depth, d, whilst an acoustic time-scale, dj (i?*T ) ' (where 
i?* is the gas constant) is used to rescale time. After rescaling p, T, u and B in an appropriate way (see 
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Parameter 


Definition 


Values used 


7 


Cp/c V 


5/3 


m 


gd/R*AT - 1 


1.0 


9 


at/Tq 


10.0 


K 


K/podcpiR^Ta) 1 ' 2 


0.2 


Co 


TjpoCp/K 


0.1 


a 


licpjK 


1.0 


F 


B 2 /po[j.qR*T 


Variable 



Table 1: The non-dimensional parameters in the governing equations for compressible magnetoconvection. 
Note that po corresponds to the unperturbed density at the upper surface of the layer. All other parameters 
are as defined in the text. 



Bushbv fc Houghton! (|2005f ) for more details), the governing equations for this system can be written in the 
following form: 



V ■ (pu) , (1) 



dp 
dt ' 

^ (pu) = -V( P T + F\B\ 2 /2) + 9(m + l)pz 
+ V • (F BB - puu + hot) , 

V x (u x B - kCqV x B) , (3) 



<9B 

~8t 



^ = - u • VT - ( 7 - 1) TV • u + ?1 V 2 T 
at p 



k(7-1) 



P 

where the components of the stress-tensor, t, are given by 



(4) 

a^^ + FColVxBI 2 ), 



_ duj | duj 2 du k x . . 

Tij ~d Xj + d Xi 3dx k ( ) 

whilst B satisfies the standard constraint that V • B = 0. The pressure P satisfies 

P = P T. (6) 

Equation (1) describes the conservation of mass, whilst Equation (2) is simply the momentum equation, written 
in conservative form. Note that the 6(m + l)pz term in Equation (2) represents the effects of gravity whilst 
the two terms that are quadratic in B correspond to the Lorentz force. The evolution of the magnetic field 
is governed by the standard magnetic induction equation (3). The final two terms in the thermal equation 
(4) represent the effects of viscous and ohmic heating. The non-dimensional parameters that appear in the 
governing equations are defined in Table 1. For later convenience, we also introduce the Chandrasekhar number, 

Q = F/k 2 (o<t, (7) 

and the Rayleigh number, 

Ra = (m + 1 - m 7 )(l + 8/2) 2 " 1 - 1 (to ± 1)0 . (8) 

The Chandrasekhar number is a measure of the strength of the imposed magnetic field whilst the Rayleigh 
number measures the destabilising influence of the temperature gradient that is imposed across the layer. If 
all other parameters are fixed, varying Ra and Q is equivalent to varying F and k. 

There is a non-trivial equilibrium solution to these governing equations, corresponding to a static polytropic 
layer with a uniform, vertical magnetic field: 

u = 0, T(z) = l + 6z, p{z) = (1 + 6z) m , B = 2. (9) 
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Figure 1: The temperature distribution in a horizontal plane just below the upper surface of the computational 
domain for Q = 120 (top) and Q = 150 (bottom). Brighter contours correspond to regions of warmer fluid. 
These simulations were both evolved from a static state. 



In this equilibrium state, the parameters that are given in Table 1 imply that the layer of fluid is highly 
stratified, with the temperature and density both varying by a factor of 11 across the depth of the domain. 
When there is no magnetic field present (i.e. when Q — or, equivalently, F = 0), the critical Rayleigh 
number for the onset of convection in this layer is approximately, Ra = 1189. For the parameters given in 
Table 1, Ra = 6000, hence this lay er is convectively unst able in the absence of an imposed magnetic field. As 
in Boussinesq magnetoconvection ( Chandrasekhar . 196lh . magnetic fields tend to inhibit convection. Hence, 
non-zero values of Q lead to an increase in the critical Rayleigh number. Table 1 also gives values for two 
diffusivity ratios, a and For simplicity, we set a — 1.0. The value of £o plays a crucial role in determining the 
near-onset bifurcation structure. As was the case in the truncated Boussinesq model of Blanch-flower ( 19991) . a 
choice of £o = 0.1 ensures that the equilibrium state can be unstable to oscillatory modes of convection as well 
as stationary modes. 



3 Numerical results 

Given the complexity of the governing equations, it is necessary to solve these numerically. We cover the 
Cartesian domain with a computational mesh, typically consisting of 128 x 128 x 96 grid points. Using 
standard Fast Fourier Transform (FFT) libraries, all horizontal derivatives are evaluated in Fourier space. 
Fourth-order finite differences are used to calculate the vertical derivatives. The temporal evolution of this 
system is determined by an explicit third-order Adams-Bashforth scheme. This code is parallelised using MPI. 



3.1 Supercritical convection 

As described in Section [2j one simple solution of these governing equations corresponds to a static polytropic 
layer in the presence of a uniform, vertical magnetic field (equation [5]). We use this static solution (plus a small 
thermal perturbation) as an initial condition for the code. For the parameter values that are given in Table 1, 
this initial state is convectively unstable provided that the Chandrasekhar number does not exceed a value of 
(approximately) Q = 160. Since this equilibrium solution is unstable when Q < 160, we refer to this region of 
parameter space as the 'supercritical' parameter regime. 

Figure Q] shows snapshots of the system for two different values of Q (Q = 120 and Q = 150), once a 
statistically-steady state has been reached. Each plot shows the temperature distribution across a horizontal 
layer, just below the upper surface of the domain. When Q = 120, the solution is characterised by several 
time-dependent plumes from which most of the surface magnetic field has been expelled by the convective 
motions. Elsewhere the flow is dominated by magnetic effects, so much so that the Lorentz force is strong 
enough to (almost) comple tely inhibit convect ion. Solutions of this form are reminiscent of the 'flux-separated' 
states that were found by IWeiss et al. (l2002h although, in that study, small-scale convective cells tended to 
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a 2 * t b 

Figure 2: Like the lower part of Figure Q] this shows the temperature distribution, for Q = 150, across a 
horizontal layer just below the upper surface of the domain. This solution has been generated by following the 
flux-separated solution branch. 



be found in the magnetically-dominated regions. Here the magnetic suppression of convection is much more 
pronounced in this flux-separated state. A very different form of solution is found when Q = 150 (lower part 
of Fig. [1} . In this simulation, the convection forms an oscillatory pattern that is distributed across the whole 

of the domain. _ ^ 

Linear theory (see, for example, IChandrasekhaii 1 1 9 6 lh predicts oscillatory convection in the magnetically- 
dominated regime, at least for low values of Co- Hence the oscillatory behaviour of the simulation at Q = 150 
is unsurprising. The convection at Q = 120 is also time-dependent, as noted above, but this flux-separated 
state is representative of a different solution branch. Although th e present model is m uch more complicated 
than the two-dimensional Boussinesq model that was considered bv iBlanchflower ( 1999t h it is clearly of interest 
to relate the two studies. The flux separated state at Q = 12 would correspond to one of the multiple roll 
states that are shown in the bifurcation diagram in Figure 3 of lBlanchflower ( 19991 ). The existence of different 
solution branches raises the possibility that a given set of para meter values could be associated with more than 
one stable state. Following the procedure that is described bv IBlanchflower ( 19991) . solution branches can be 
followed by starting a simulation with fully developed convection, rather than evolving it from a static state. 
Having adjusted the value of Q, the simulation can then be evolved again until a statistically-steady state has 
been found. For example, if the Q = 120 solution is taken as an initial condition, the flux-separated solution 
branch can be followed by gradually increasing the strength of the the magnetic field (by increasing the value 
of Q). The outcome of such a procedure is illustrated in Figure [2] which shows the temperature distribution 
for a flux-separated solution at Q — 150. Comparing this plot with the lower part of Figure [TJ it is clear that 
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(at least) two distinct solutions exist for this set of parameter values. This flux-separated solution has smaller 
convective plumes than the Q = 120 case. This is a consequence of the fact that magnetic fields are now 
strong enough to inhibit convection in a larger proportion of the computational domain. The typical scale of 
convection is reduced by increasing the strength of the imposed magnetic field. Additionally the rate at which 
convective plumes merge together, and split apart, reduces as the strength of the imposed magnetic field is 
increased. 



3.2 The subcritical parameter regime 

As noted above, the static polytropic layer is (linearly) stable to convective perturbations for values of the 
Chandrasekhar number in excess of approximately Q = 160. Hence, we refer to this range of parameter space 
as the 'subcritical' regime. In order to find non-trivial behaviour in this parameter regime, it is clearly necessary 
to adopt a non-static initial condition for any simulations that are carried out. Given the results that were 
presented in the previous section, it is natural to try to track the flux-separated solution branch into this 
subcritical regime by incrementally increasing the value of Q ( as before, look i ng for a statistically-steady state 
before each increment). By following an analogous procedure. iBlanchflowerl (|l999l ) found localised states in a 
two-dimensional model, so this would appear to be the most sensible approach. 

As the value of Q is gradually increased, following the flux-separated branch, the dynamical influence of 
the magnetic field becomes greater. This leads to a reduction in both the number and scale of the field-free 
convective plumes. This process continues until a single steady plume remains. Thi s plume is the three - 
dimensional analogue of the two-dimensional 'convecton' solutions that were found bv iBlanchflower ( 1999t ). 



This three-dimensional convecton is illustrated in Figures and |U At Q = 215, this localised convective 
plume is almost axisymmetric, being slightly elongated in the x-direction, with a broad central upwelling, 
surrounded by a narrow downflow region. Convection is completely suppressed everywhere else by a strong, 
uniform, vertical magnetic field. The magnetic field distribution within the convecton is more complicated. 
At the surface, the field is almost completely expelled by the diverging convective flows. Towards the base of 
the plume, converging convective flows lead to an accumulation of vertical magnetic flux at the base of the 
convective upflow. The stratification of the layer clearly does play a role in determining the structure of this 
localised convective feature: slightly larger temperature perturbations are found near the top of the layer. It 
should be stressed that this convecton is steady and is therefore not simply a transient phenomenon. It is also 
worth noting that, although the horizontally averaged convective flux is small, the local perturbations to the 
thermodynamic variables are of order unity within the convecton itself. Therefore this is also a dynamically 
significant feature, albeit a highly localised one. 

It should be noted that this convecton is not truly localised rather, by virtue of the horizontal boundary 
conditions, it is part of a periodic array of such convective structures. As seen in figure U] each structure is 
separated by approximately 6 non-dimensional units which is approximately 3 times the width of the convective 
structure. Therefore, any influence from neighbouring convectons will be small. An initial investigation in larger 
computational domains confirms this to be the case. Similar localised states can be found, although the range 
of Q for which they exist does change with domain size. Attempts to find an analytic description of a perfectly 
axisymmetric convecton are ongoing. 

Having found this convecton, this solution branch can also be tracked to determine its range of stability. 
In fact, the solution that is shown in Figures [3] and [4] corresponds to the largest value of Q at which the 
convecton is found to be stable. This is also the value of Q at which the greatest degree of localisation is 
found. As this solution branch is followed into the weaker field regime (by gradually decreasing the value of 
Q), the convecton grows, becoming increasingly asymmetric. For example, at Q = 190 the solution is still 
highly localised and steady, but the convective cell is about 25% wider than the convecton at Q = 215 and 
is more elongated in the x-direction. For this larger localised state, it is plausible that the finite size of the 
computational domain is (weakly) influencing the symmetry of the convective plume. Reducing Q still further, 
we find that the solution is still mostly localised, although there are now weak (but significant) fluctuations 
elsewhere in the domain. More importantly, the convecton is no longer steady. Instead, the plume 'wobbles' 
periodically about its central axis (as illustrated in Figure [SJ. This is not a true oscillatory state, and the 
amplitude of the fluctuations is small, but this behaviour is certainly interesting. Sim ilar vacillation has bee n 
seen as a way for steady convection patterns to lose stability via a Hopf bifurcation ( Rucklidge et al. . 19931) . 



It is also possible that the observed oscillations are a consequence of the finite geometry, in the sense that the 
convecton could be interacting with periodic copies of itself. Calculations in larger domains to investigate this 
further are currently underway. 



G 



a 2 + £ B 

Figure 3: Steady, localised convection at Q=215. This shows the temperature distribution in a horizontal 
plane just below the upper surface of the computational domain. Taking into account the periodic boundary 
conditions, this solution corresponds to a single convective plume. The convecton is close to axisymmetric, but 
slightly elongated in the x-direction. 
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Figure 4: As Figure El these plots show snapshots of a convecton at Q = 215. Top: The perturbation to the 
background temperature distribution along a vertical slice through x — 7.5. Contours range from —0.11 (black) 
to 0.87 (white). Bottom: The vertical component of the magnetic held in the same vertical slice. Contours 
range from 0.22 (black) to 3.67 (red). 
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Figure 5: Snapshots of convection along the convecton solution branch (at Q — 160). These plots show the 
temperature distribution in a horizontal plane just below the upper surface of the computational domain. 
These plots are separated by 4.55 time units, the full cycle period is 18.20 time units. 
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4 Discussion and Conclusions 



In these numerical simulations we have found strongly localised convective plumes in the subcritical regime 
of compressible magnetoconvection. These states are not embedded within a background of diffuse weak 
convection, but rather a static fluid (see figure It is believed that strongly localised, steady states of this 
type have not previously been observed in three spatial dimensions. Whether umbral dots ever exhibit such an 
extreme degree of localisation is unclear, but these results do suggest that highly localised plumes could occur 
if the magnetic fields within sunspot umbrae are strong enough (locally) to produce subcritical conditions for 
magnetoconvection. 

The strongly localised states found in this work ar e found by continu ation of flux-separated states to higher 
values of Q. A similar procedure was carried out by IWeiss et all (|2002l ). however, in their work it was found 
that for sufficiently strong impo s ed ma gnetic fields the flux separated solution was lost in favour of small-scale 
regular convection. I Weiss et al.l (|2002f) used the parameter value Co = 0.2 whereas in this work Co = 0.1. It is 
well-know n that the value of C n controls the near-onset bifurcation structure. Comparison of this work with the 
findings of IWeiss et al. ( 2002 ) indicate that it also controls whether or not the solution branch corresponding 
to flux-separated states continues into the subcritical regime. 

Although it is appealing to relate these highly localised convective plumes to umbral dots, one aspect of 
these solutions that does not compare favourably with umbral dots is that fact that these localised states are 
steady (or, at best, weakly time-dependent, as shown in Figure [5]). So, unlike the observed umbral dots, these 
solutions do not have finite lifetimes. This may be associated with the simplified nature of the model, which 
assumes a uniform background state: spatiotemporal inhomogeneities in the umbral background may play an 
important role in determining the lifetime of an umbral dot. Alternatively, time-dependent behaviour may 
be the result of interacting plumes. Perhaps the time-dependent state that is shown in Figure [2] is a more 
realistic representation of the distribution of convective plumes within sunspot umbrae. Work is ongoing to 
establish whether or no t it is possible to find truly localised states that exhibit significant time-dependence. 
The simplified model of Blanchflower fc Weiss! ( 2002 ) suggests that oscillatory solutions should exist, however, 
while we cannot rule out this possibility, oscillatory localised states have not yet been observed in the full 
system of equations. 

This model is clearly a highly idealised representation of magnetoconvection within sunspot umbrae. Truly 
localised states have not yet been found in more realistic models of photospheric magnetoconvection. It would 
be of great interest to establish whether or not processes such as radiative transfer promote (or inhibit) the 
formation of strongly localised plumes. It should also be stressed that the range of parameter values that can be 
considered in numerical models of this type bear little resemblance to true photospheric values (although this 
is true of all numerical magnetoconvection calculations, regardless of the level of physical complexity). Having 
said that, we believe that results from models of this type do provide some useful insights into photospheric 
magnetoconvection. 

Magnetoconvection is not the only branch of physics in which strongly localised state s are relevant. Similar 
state s have now been found in a variety of other s ystems, including binary fluid convection ( Bergeon fc KnoblochL 
2008), buckling rod pro blems, nonlinear optics (Vladimirov et al. . 20021 ) and experiments on a ferrofluid in an 
applied magnetic field (jRichter fc Barashenkovl 120051 ) . In a more ab stract setting, localised states h ave also 
been found in the one-dimensional bistable Swift-Hohenberg equation ( Burke fc Knoblo"chl ,2006, 200730) • The 
results presented in this paper share many similarities with the behaviour of the Swift-Hohenberg equation 
in one dimension, as well as the other physical systems mentioned above. There is now a good theoretical 
understanding of localised states in one extended dimension, howev er, in higher dime nsions the problem is not 
well understood, and has recently been posed as an open proble m ( Knoblocb . 20081) . Having said that, some 
recent progress has been made. For example, iLlovd et al. (2008) have investigated the properties of spatially 
localised states in the two-dimensional Swift-Hohenberg equation. In future work, we intend to further ex- 
plore the relationship between the behaviour of this simple pattern forming system and the magnetoconvection 
equations. 
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